{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15\n"
     ]
    }
   ],
   "source": [
    "import sys,os\n",
    "add_path = os.path.abspath('../python')\n",
    "if not add_path in sys.path:\n",
    "    sys.path.append(add_path)\n",
    "print(sys.path.index(add_path))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "from numpy import vstack, hstack\n",
    "from segbench import csaAssign, correspondPixels, matchEdgeMaps2D, matchEdgeMaps3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  5.  1.]\n",
      " [ 2.  6.  1.]\n",
      " [ 3.  4.  5.]]\n"
     ]
    }
   ],
   "source": [
    "data = np.array([[1, 4, 3],\n",
    "    [2, 5, 3],\n",
    "    [3, 6, 3],\n",
    "    [1, 5, 1],\n",
    "    [2, 6, 1],\n",
    "    [3, 4, 5]])\n",
    "n = 6\n",
    "e1 = csaAssign(n,data)\n",
    "print e1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.273e-02s\n",
      "1000 500250 [  4.00000000e+00   1.00500000e+03   1.00000000e+00]\n"
     ]
    }
   ],
   "source": [
    "# Unit test for csaAssign.\n",
    "n = 1000\n",
    "dg = np.random.rand(n,n)>.5\n",
    "m = np.sum(dg.flat)\n",
    "i = np.array(dg.ravel().nonzero()).astype(float)\n",
    "g1 = 1+np.floor(np.divide(i,n))\n",
    "g2 = (1+np.floor(np.divide(i,n)),\n",
    "     1+ np.mod(i,n)+n,\n",
    "     1+np.floor(np.random.rand(1,m))*1000)\n",
    "g2 = np.array(np.vstack(g2)).T\n",
    "t = time.time()\n",
    "e2 = csaAssign(n*2,g2)\n",
    "print '{:.3e}s'.format(time.time() - t)\n",
    "esum = np.sum(e2,axis=0)\n",
    "if not all(esum[0:2]==(n*(n+1)/2,n*(n+1)/2+n*n)): raise ValueError(\"bug\")\n",
    "if not np.sum(esum[0:2])==(2*n*(2*n+1)/2): raise ValueError(\"bug\")\n",
    "print n,m, e2[3][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n=1 thresh=0.00e+00 4.29e-05s 2 2 0\n",
      "n=1 thresh=1.11e-01 3.60e-05s 2 2 0\n",
      "n=1 thresh=2.22e-01 2.38e-05s 2 2 0\n",
      "n=1 thresh=3.33e-01 2.29e-05s 2 2 0\n",
      "n=1 thresh=4.44e-01 1.22e-05s 1 0 1\n",
      "n=1 thresh=5.56e-01 2.81e-05s 1 0 1\n",
      "n=1 thresh=6.67e-01 1.31e-05s 1 0 1\n",
      "n=1 thresh=7.78e-01 1.22e-05s 0 0 0\n",
      "n=1 thresh=8.89e-01 1.00e-05s 0 0 0\n",
      "n=1 thresh=1.00e+00 9.06e-06s 0 0 0\n",
      "n=2 thresh=0.00e+00 4.79e-05s 8 8 0\n",
      "n=2 thresh=1.11e-01 4.51e-05s 8 8 0\n",
      "n=2 thresh=2.22e-01 4.89e-05s 8 8 0\n",
      "n=2 thresh=3.33e-01 2.60e-05s 4 4 0\n",
      "n=2 thresh=4.44e-01 2.50e-05s 6 4 2\n",
      "n=2 thresh=5.56e-01 2.10e-05s 4 2 2\n",
      "n=2 thresh=6.67e-01 1.00e-05s 2 0 2\n",
      "n=2 thresh=7.78e-01 1.00e-05s 0 0 0\n",
      "n=2 thresh=8.89e-01 1.00e-05s 2 0 2\n",
      "n=2 thresh=1.00e+00 1.00e-05s 0 0 0\n",
      "n=3 thresh=0.00e+00 1.10e-04s 18 18 0\n",
      "n=3 thresh=1.11e-01 1.05e-04s 18 18 0\n",
      "n=3 thresh=2.22e-01 7.20e-05s 15 12 3\n",
      "n=3 thresh=3.33e-01 4.39e-05s 10 6 4\n",
      "n=3 thresh=4.44e-01 2.69e-05s 10 4 6\n",
      "n=3 thresh=5.56e-01 2.60e-05s 9 4 5\n",
      "n=3 thresh=6.67e-01 1.00e-05s 3 0 3\n",
      "n=3 thresh=7.78e-01 2.19e-05s 6 2 4\n",
      "n=3 thresh=8.89e-01 1.10e-05s 1 0 1\n",
      "n=3 thresh=1.00e+00 1.10e-05s 0 0 0\n",
      "n=5 thresh=0.00e+00 3.58e-04s 50 50 0\n",
      "n=5 thresh=1.11e-01 2.77e-04s 46 42 4\n",
      "n=5 thresh=2.22e-01 1.80e-04s 37 24 13\n",
      "n=5 thresh=3.33e-01 1.59e-04s 32 22 10\n",
      "n=5 thresh=4.44e-01 1.26e-04s 29 18 11\n",
      "n=5 thresh=5.56e-01 8.08e-05s 23 12 11\n",
      "n=5 thresh=6.67e-01 4.01e-05s 16 6 10\n",
      "n=5 thresh=7.78e-01 1.38e-05s 8 0 8\n",
      "n=5 thresh=8.89e-01 1.19e-05s 7 0 7\n",
      "n=5 thresh=1.00e+00 2.81e-05s 0 0 0\n",
      "n=8 thresh=0.00e+00 1.16e-03s 128 128 0\n",
      "n=8 thresh=1.11e-01 8.69e-04s 112 102 10\n",
      "n=8 thresh=2.22e-01 7.25e-04s 99 78 21\n",
      "n=8 thresh=3.33e-01 6.64e-04s 99 76 23\n",
      "n=8 thresh=4.44e-01 1.72e-04s 61 28 33\n",
      "n=8 thresh=5.56e-01 1.70e-04s 54 24 30\n",
      "n=8 thresh=6.67e-01 9.89e-05s 40 14 26\n",
      "n=8 thresh=7.78e-01 6.70e-05s 27 10 17\n",
      "n=8 thresh=8.89e-01 2.81e-05s 16 2 14\n",
      "n=8 thresh=1.00e+00 1.50e-05s 0 0 0\n",
      "n=13 thresh=0.00e+00 3.08e-03s 338 338 0\n",
      "n=13 thresh=1.11e-01 2.59e-03s 300 268 32\n",
      "n=13 thresh=2.22e-01 1.66e-03s 249 178 71\n",
      "n=13 thresh=3.33e-01 1.07e-03s 208 126 82\n",
      "n=13 thresh=4.44e-01 1.07e-03s 197 106 91\n",
      "n=13 thresh=5.56e-01 6.22e-04s 156 80 76\n",
      "n=13 thresh=6.67e-01 4.66e-04s 127 52 75\n",
      "n=13 thresh=7.78e-01 8.32e-05s 63 10 53\n",
      "n=13 thresh=8.89e-01 4.29e-05s 27 4 23\n",
      "n=13 thresh=1.00e+00 2.19e-05s 0 0 0\n",
      "n=22 thresh=0.00e+00 1.32e-02s 968 966 2\n",
      "n=22 thresh=1.11e-01 1.01e-02s 871 790 81\n",
      "n=22 thresh=2.22e-01 8.18e-03s 752 596 156\n",
      "n=22 thresh=3.33e-01 5.05e-03s 641 430 211\n",
      "n=22 thresh=4.44e-01 3.61e-03s 534 304 230\n",
      "n=22 thresh=5.56e-01 1.82e-03s 403 168 235\n",
      "n=22 thresh=6.67e-01 9.48e-04s 318 116 202\n",
      "n=22 thresh=7.78e-01 4.60e-04s 217 54 163\n",
      "n=22 thresh=8.89e-01 4.51e-05s 106 0 106\n",
      "n=22 thresh=1.00e+00 4.29e-05s 0 0 0\n",
      "n=36 thresh=0.00e+00 5.10e-02s 2592 2584 8\n",
      "n=36 thresh=1.11e-01 3.52e-02s 2287 2002 285\n",
      "n=36 thresh=2.22e-01 2.60e-02s 2034 1584 450\n",
      "n=36 thresh=3.33e-01 1.56e-02s 1693 1080 613\n",
      "n=36 thresh=4.44e-01 1.15e-02s 1439 790 649\n",
      "n=36 thresh=5.56e-01 7.29e-03s 1157 494 663\n",
      "n=36 thresh=6.67e-01 3.44e-03s 871 298 573\n",
      "n=36 thresh=7.78e-01 1.25e-03s 598 130 468\n",
      "n=36 thresh=8.89e-01 3.11e-04s 288 30 258\n",
      "n=36 thresh=1.00e+00 9.61e-05s 0 0 0\n",
      "n=60 thresh=0.00e+00 1.56e-01s 7200 7182 18\n",
      "n=60 thresh=1.11e-01 1.24e-01s 6420 5702 718\n",
      "n=60 thresh=2.22e-01 8.80e-02s 5633 4386 1247\n",
      "n=60 thresh=3.33e-01 5.59e-02s 4798 3176 1622\n",
      "n=60 thresh=4.44e-01 4.04e-02s 4088 2316 1772\n",
      "n=60 thresh=5.56e-01 2.00e-02s 3137 1408 1729\n",
      "n=60 thresh=6.67e-01 1.02e-02s 2337 790 1547\n",
      "n=60 thresh=7.78e-01 5.29e-03s 1653 374 1279\n",
      "n=60 thresh=8.89e-01 8.93e-04s 779 76 703\n",
      "n=60 thresh=1.00e+00 2.34e-04s 0 0 0\n",
      "n=100 thresh=0.00e+00 6.19e-01s 20000 19946 54\n",
      "n=100 thresh=1.11e-01 4.29e-01s 17750 15718 2032\n",
      "n=100 thresh=2.22e-01 3.36e-01s 15588 12110 3478\n",
      "n=100 thresh=3.33e-01 1.99e-01s 13332 8848 4484\n",
      "n=100 thresh=4.44e-01 1.49e-01s 11118 6206 4912\n",
      "n=100 thresh=5.56e-01 7.63e-02s 8844 3902 4942\n",
      "n=100 thresh=6.67e-01 4.24e-02s 6633 2176 4457\n",
      "n=100 thresh=7.78e-01 1.61e-02s 4513 1072 3441\n",
      "n=100 thresh=8.89e-01 3.20e-03s 2242 240 2002\n",
      "n=100 thresh=1.00e+00 6.41e-04s 0 0 0\n"
     ]
    }
   ],
   "source": [
    "# Unit test for correspondPixels.\n",
    "for n in np.round(np.logspace(0,2,10)).astype(int):\n",
    "    for thresh in np.linspace(0,1,10):\n",
    "        print 'n={:d} thresh={:1.2e}'.format(n,thresh),\n",
    "        bmap1 = (np.random.rand(n,n)>thresh).astype(float)\n",
    "        bmap2 = (np.random.rand(n,n)>thresh).astype(float)\n",
    "        t = time.time()\n",
    "        cost,oc,match1,match2 = correspondPixels(bmap1,bmap2,.005)\n",
    "        print '{:1.2e}s'.format(time.time() - t),\n",
    "        npix = len(bmap1.nonzero()[0]) + len(bmap2.nonzero()[0])\n",
    "        ncor = len(match1.nonzero()[0]) + len(match1.nonzero()[0])\n",
    "        nout = len(((bmap1>0)&(match1==0)).nonzero()[0]) + len(((bmap2>0)&(match2==0)).nonzero()[0])\n",
    "        print npix,ncor,nout\n",
    "        if ncor+nout!=npix: raise ValueError(\"bug\")\n",
    "\n",
    "        match1f = match1.ravel()\n",
    "        idx1 = match1f.astype(int)\n",
    "        a1 = match1f.nonzero()[0]\n",
    "        a2 = match1f[idx1.nonzero()[0]].ravel()\n",
    "        a = vstack((a1,a2)).T\n",
    "\n",
    "        match2f = match2.ravel()\n",
    "        idx2 = match2f.astype(int)\n",
    "        b1 = match2f.nonzero()[0]\n",
    "        b2 = match2f[idx2.nonzero()[0]].ravel()\n",
    "        b = vstack((b2,b1)).T\n",
    "        b = b[b[:,0].argsort(),]\n",
    "        if np.sum(a-b): raise ValueError(\"bug\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n=1 thresh=0.00e+00 4.60e-05s 2 + 0 = 2 == 2\n",
      "n=1 thresh=1.11e-01 1.60e-05s 0 + 1 = 1 == 1\n",
      "n=1 thresh=2.22e-01 1.91e-05s 0 + 1 = 1 == 1\n",
      "n=1 thresh=3.33e-01 3.91e-05s 2 + 0 = 2 == 2\n",
      "n=1 thresh=4.44e-01 2.10e-05s 0 + 1 = 1 == 1\n",
      "n=1 thresh=5.56e-01 2.10e-05s 0 + 0 = 0 == 0\n",
      "n=1 thresh=6.67e-01 2.19e-05s 0 + 1 = 1 == 1\n",
      "n=1 thresh=7.78e-01 1.91e-05s 0 + 1 = 1 == 1\n",
      "n=1 thresh=8.89e-01 1.98e-05s 0 + 0 = 0 == 0\n",
      "n=1 thresh=1.00e+00 2.00e-05s 0 + 0 = 0 == 0\n",
      "n=2 thresh=0.00e+00 6.60e-05s 8 + 0 = 8 == 8\n",
      "n=2 thresh=1.11e-01 3.10e-05s 4 + 1 = 5 == 5\n",
      "n=2 thresh=2.22e-01 3.70e-05s 6 + 0 = 6 == 6\n",
      "n=2 thresh=3.33e-01 3.00e-05s 4 + 2 = 6 == 6\n",
      "n=2 thresh=4.44e-01 1.10e-05s 0 + 4 = 4 == 4\n",
      "n=2 thresh=5.56e-01 3.79e-05s 2 + 1 = 3 == 3\n",
      "n=2 thresh=6.67e-01 3.00e-05s 4 + 1 = 5 == 5\n",
      "n=2 thresh=7.78e-01 1.69e-05s 0 + 0 = 0 == 0\n",
      "n=2 thresh=8.89e-01 2.31e-05s 0 + 2 = 2 == 2\n",
      "n=2 thresh=1.00e+00 2.79e-05s 0 + 0 = 0 == 0\n",
      "n=3 thresh=0.00e+00 1.90e-04s 18 + 0 = 18 == 18\n",
      "n=3 thresh=1.11e-01 1.41e-04s 18 + 0 = 18 == 18\n",
      "n=3 thresh=2.22e-01 6.29e-05s 8 + 5 = 13 == 13\n",
      "n=3 thresh=3.33e-01 6.91e-05s 8 + 4 = 12 == 12\n",
      "n=3 thresh=4.44e-01 3.91e-05s 6 + 5 = 11 == 11\n",
      "n=3 thresh=5.56e-01 4.70e-05s 8 + 3 = 11 == 11\n",
      "n=3 thresh=6.67e-01 2.31e-05s 2 + 5 = 7 == 7\n",
      "n=3 thresh=7.78e-01 1.00e-05s 0 + 1 = 1 == 1\n",
      "n=3 thresh=8.89e-01 1.00e-05s 0 + 0 = 0 == 0\n",
      "n=3 thresh=1.00e+00 9.06e-06s 0 + 0 = 0 == 0\n",
      "n=5 thresh=0.00e+00 4.19e-04s 50 + 0 = 50 == 50\n",
      "n=5 thresh=1.11e-01 2.74e-04s 38 + 6 = 44 == 44\n",
      "n=5 thresh=2.22e-01 1.91e-04s 28 + 9 = 37 == 37\n",
      "n=5 thresh=3.33e-01 1.42e-04s 18 + 12 = 30 == 30\n",
      "n=5 thresh=4.44e-01 9.49e-05s 12 + 14 = 26 == 26\n",
      "n=5 thresh=5.56e-01 3.10e-05s 4 + 12 = 16 == 16\n",
      "n=5 thresh=6.67e-01 3.81e-05s 6 + 9 = 15 == 15\n",
      "n=5 thresh=7.78e-01 1.29e-05s 0 + 9 = 9 == 9\n",
      "n=5 thresh=8.89e-01 1.10e-05s 0 + 6 = 6 == 6\n",
      "n=5 thresh=1.00e+00 1.10e-05s 0 + 0 = 0 == 0\n",
      "n=8 thresh=0.00e+00 1.16e-03s 128 + 0 = 128 == 128\n",
      "n=8 thresh=1.11e-01 8.21e-04s 94 + 17 = 111 == 111\n",
      "n=8 thresh=2.22e-01 5.93e-04s 74 + 22 = 96 == 96\n",
      "n=8 thresh=3.33e-01 4.86e-04s 66 + 26 = 92 == 92\n",
      "n=8 thresh=4.44e-01 2.42e-04s 34 + 30 = 64 == 64\n",
      "n=8 thresh=5.56e-01 1.79e-04s 26 + 32 = 58 == 58\n",
      "n=8 thresh=6.67e-01 1.00e-04s 12 + 29 = 41 == 41\n",
      "n=8 thresh=7.78e-01 3.48e-05s 4 + 31 = 35 == 35\n",
      "n=8 thresh=8.89e-01 2.88e-05s 2 + 9 = 11 == 11\n",
      "n=8 thresh=1.00e+00 1.41e-05s 0 + 0 = 0 == 0\n",
      "n=13 thresh=0.00e+00 3.90e-03s 338 + 0 = 338 == 338\n",
      "n=13 thresh=1.11e-01 2.73e-03s 264 + 36 = 300 == 300\n",
      "n=13 thresh=2.22e-01 1.72e-03s 202 + 60 = 262 == 262\n",
      "n=13 thresh=3.33e-01 1.38e-03s 150 + 70 = 220 == 220\n",
      "n=13 thresh=4.44e-01 8.75e-04s 98 + 90 = 188 == 188\n",
      "n=13 thresh=5.56e-01 6.17e-04s 72 + 79 = 151 == 151\n",
      "n=13 thresh=6.67e-01 3.30e-04s 48 + 77 = 125 == 125\n",
      "n=13 thresh=7.78e-01 1.23e-04s 14 + 61 = 75 == 75\n",
      "n=13 thresh=8.89e-01 4.91e-05s 6 + 33 = 39 == 39\n",
      "n=13 thresh=1.00e+00 2.19e-05s 0 + 0 = 0 == 0\n",
      "n=22 thresh=0.00e+00 1.29e-02s 966 + 2 = 968 == 968\n",
      "n=22 thresh=1.11e-01 8.75e-03s 776 + 92 = 868 == 868\n",
      "n=22 thresh=2.22e-01 6.87e-03s 572 + 178 = 750 == 750\n",
      "n=22 thresh=3.33e-01 7.69e-03s 446 + 215 = 661 == 661\n",
      "n=22 thresh=4.44e-01 3.40e-03s 300 + 240 = 540 == 540\n",
      "n=22 thresh=5.56e-01 1.62e-03s 174 + 241 = 415 == 415\n",
      "n=22 thresh=6.67e-01 9.91e-04s 100 + 224 = 324 == 324\n",
      "n=22 thresh=7.78e-01 3.65e-04s 40 + 170 = 210 == 210\n",
      "n=22 thresh=8.89e-01 1.41e-04s 12 + 88 = 100 == 100\n",
      "n=22 thresh=1.00e+00 4.32e-05s 0 + 0 = 0 == 0\n",
      "n=36 thresh=0.00e+00 4.63e-02s 2590 + 2 = 2592 == 2592\n",
      "n=36 thresh=1.11e-01 3.19e-02s 2044 + 265 = 2309 == 2309\n",
      "n=36 thresh=2.22e-01 2.85e-02s 1560 + 465 = 2025 == 2025\n",
      "n=36 thresh=3.33e-01 1.79e-02s 1190 + 584 = 1774 == 1774\n",
      "n=36 thresh=4.44e-01 1.19e-02s 830 + 641 = 1471 == 1471\n",
      "n=36 thresh=5.56e-01 5.35e-03s 484 + 624 = 1108 == 1108\n",
      "n=36 thresh=6.67e-01 3.15e-03s 280 + 595 = 875 == 875\n",
      "n=36 thresh=7.78e-01 1.66e-03s 172 + 458 = 630 == 630\n",
      "n=36 thresh=8.89e-01 3.42e-04s 30 + 230 = 260 == 260\n",
      "n=36 thresh=1.00e+00 9.61e-05s 0 + 0 = 0 == 0\n",
      "n=60 thresh=0.00e+00 1.51e-01s 7180 + 20 = 7200 == 7200\n",
      "n=60 thresh=1.11e-01 1.15e-01s 5662 + 732 = 6394 == 6394\n",
      "n=60 thresh=2.22e-01 7.30e-02s 4278 + 1296 = 5574 == 5574\n",
      "n=60 thresh=3.33e-01 5.78e-02s 3198 + 1619 = 4817 == 4817\n",
      "n=60 thresh=4.44e-01 3.89e-02s 2292 + 1746 = 4038 == 4038\n",
      "n=60 thresh=5.56e-01 2.13e-02s 1460 + 1707 = 3167 == 3167\n",
      "n=60 thresh=6.67e-01 1.08e-02s 816 + 1595 = 2411 == 2411\n",
      "n=60 thresh=7.78e-01 4.02e-03s 340 + 1192 = 1532 == 1532\n",
      "n=60 thresh=8.89e-01 1.28e-03s 90 + 698 = 788 == 788\n",
      "n=60 thresh=1.00e+00 2.18e-04s 0 + 0 = 0 == 0\n",
      "n=100 thresh=0.00e+00 5.91e-01s 19968 + 32 = 20000 == 20000\n",
      "n=100 thresh=1.11e-01 3.61e-01s 15846 + 1982 = 17828 == 17828\n",
      "n=100 thresh=2.22e-01 2.71e-01s 12058 + 3497 = 15555 == 15555\n",
      "n=100 thresh=3.33e-01 1.95e-01s 8946 + 4425 = 13371 == 13371\n",
      "n=100 thresh=4.44e-01 1.44e-01s 6164 + 4926 = 11090 == 11090\n",
      "n=100 thresh=5.56e-01 7.53e-02s 3930 + 5045 = 8975 == 8975\n",
      "n=100 thresh=6.67e-01 3.55e-02s 2152 + 4427 = 6579 == 6579\n",
      "n=100 thresh=7.78e-01 1.35e-02s 914 + 3445 = 4359 == 4359\n",
      "n=100 thresh=8.89e-01 3.65e-03s 258 + 2006 = 2264 == 2264\n",
      "n=100 thresh=1.00e+00 9.03e-04s 0 + 0 = 0 == 0\n"
     ]
    }
   ],
   "source": [
    "# Unit test for matchEdgeMaps2D.\n",
    "for n in np.round(np.logspace(0,2,10)).astype(int):\n",
    "    for thresh in np.linspace(0,1,10):\n",
    "        print 'n={:d} thresh={:1.2e}'.format(n,thresh),\n",
    "        bmap1 = (np.random.rand(n,n)>thresh).astype(float)\n",
    "        bmap2 = (np.random.rand(n,n)>thresh).astype(float)\n",
    "        t = time.time()\n",
    "        cost,match1,match2 = matchEdgeMaps2D(bmap1,bmap2,.005,100)\n",
    "        print '{:1.2e}s'.format(time.time() - t),\n",
    "        npix = len(bmap1.nonzero()[0]) + len(bmap2.nonzero()[0])\n",
    "        ncor = len(match1.nonzero()[0]) + len(match1.nonzero()[0])\n",
    "        nout = len(((bmap1>0)&(match1==0)).nonzero()[0]) + len(((bmap2>0)&(match2==0)).nonzero()[0])\n",
    "        print ncor,'+',nout,'=',ncor+nout,'==',npix\n",
    "        if ncor+nout!=npix: raise ValueError(\"bug\")\n",
    "        \n",
    "        match1f = match1.ravel()\n",
    "        idx1 = match1f.astype(int)\n",
    "        a1 = match1f.nonzero()[0]\n",
    "        a2 = match1f[idx1.nonzero()[0]].ravel()\n",
    "        a = vstack((a1,a2)).T\n",
    "\n",
    "        match2f = match2.ravel()\n",
    "        idx2 = match2f.astype(int)\n",
    "        b1 = match2f.nonzero()[0]\n",
    "        b2 = match2f[idx2.nonzero()[0]].ravel()\n",
    "        b = vstack((b2,b1)).T\n",
    "        b = b[b[:,0].argsort(),]\n",
    "        if np.sum(a-b): raise ValueError(\"bug\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n=1 thresh=5.00e-01 0.0000s 0 + 1 = 1 == 1\n",
      "n=1 thresh=6.33e-01 0.0000s 0 + 1 = 1 == 1\n",
      "n=1 thresh=7.67e-01 0.0000s 0 + 0 = 0 == 0\n",
      "n=1 thresh=9.00e-01 0.0000s 0 + 0 = 0 == 0\n",
      "n=4 thresh=5.00e-01 0.0003s 22 + 38 = 60 == 60\n",
      "n=4 thresh=6.33e-01 0.0002s 18 + 29 = 47 == 47\n",
      "n=4 thresh=7.67e-01 0.0000s 4 + 24 = 28 == 28\n",
      "n=4 thresh=9.00e-01 0.0000s 2 + 10 = 12 == 12\n",
      "n=16 thresh=5.00e-01 0.0308s 2014 + 2065 = 4079 == 4079\n",
      "n=16 thresh=6.33e-01 0.0176s 1098 + 1874 = 2972 == 2972\n",
      "n=16 thresh=7.67e-01 0.0054s 440 + 1443 = 1883 == 1883\n",
      "n=16 thresh=9.00e-01 0.0013s 76 + 762 = 838 == 838\n",
      "n=64 thresh=5.00e-01 7.9337s 131094 + 131138 = 262232 == 262232\n",
      "n=64 thresh=6.33e-01 2.7748s 70290 + 122430 = 192720 == 192720\n",
      "n=64 thresh=7.67e-01 0.8190s 28628 + 93645 = 122273 == 122273\n",
      "n=64 thresh=9.00e-01 0.1273s 5146 + 47371 = 52517 == 52517\n"
     ]
    }
   ],
   "source": [
    "# Unit test for matchEdgeMaps3D\n",
    "for n in np.round(np.logspace(0,6,4,base=2)).astype(int):\n",
    "    for thresh in np.linspace(.5,.9,4):\n",
    "        print 'n={:d} thresh={:1.2e}'.format(n,thresh),\n",
    "        bmap1 = (np.random.rand(n,n,n)>thresh).astype(float)\n",
    "        bmap2 = (np.random.rand(n,n,n)>thresh).astype(float)\n",
    "\n",
    "        t = time.time()\n",
    "        cost,match1,match2 = matchEdgeMaps3D(bmap1,bmap2,.000005,100)\n",
    "        print '{:1.4f}s'.format(time.time() - t),\n",
    "        npix = len(bmap1.nonzero()[0]) + len(bmap2.nonzero()[0])\n",
    "        ncor = len(match1.nonzero()[0]) + len(match1.nonzero()[0])\n",
    "        nout = len(((bmap1>0)&(match1==0)).nonzero()[0]) + len(((bmap2>0)&(match2==0)).nonzero()[0])\n",
    "        print ncor,'+',nout,'=',ncor+nout,'==',npix\n",
    "        if ncor+nout!=npix: raise ValueError(\"bug\")\n",
    "        \n",
    "        match1f = match1.ravel()\n",
    "        idx1 = match1f.astype(int)\n",
    "        a1 = match1f.nonzero()[0]\n",
    "        a2 = match1f[idx1.nonzero()[0]].ravel()\n",
    "        a = vstack((a1,a2)).T\n",
    "\n",
    "        match2f = match2.ravel()\n",
    "        idx2 = match2f.astype(int)\n",
    "        b1 = match2f.nonzero()[0]\n",
    "        b2 = match2f[idx2.nonzero()[0]].ravel()\n",
    "        b = vstack((b2,b1)).T\n",
    "        b = b[b[:,0].argsort(),]\n",
    "        if np.sum(a-b): raise ValueError(\"bug\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
